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I. INTRODUCTION 

In a non linear dispersive medium the propagation of 
a continuous wave may be impeded by the phenomena 
of Modulational Instability (MI) whereby the continuous 
wave breaks up into a train of localized pulses. This 
effect has been predicted and verified in many branches 
of physics, includingfluid dynamics 0, plasma physics F" 
non linear optics^, 0,0,0], Bose-Einstein condensates[ 

In the spectral domain the instability induces the ap- 
pearance of new frequencies whose amplitude grows ex- 
ponentially. As the dynamical process continues more 
and more energy is converted from the monochromatic 
wave to the sidebands. Simultaneously harmonics of the 
unstable frequencies appear and also grow exponentially. 

Although the process of MI is by now textbook ma- 
terial, the subsequent dynamics, and in particular the 
appearance of the higher harmonics has been relatively 
little studied. The main work we are aware about in this 
context are theoretical || and experimental investiga- 
tions in the context of the Fermi-Paska-Ulam recurrence. 
It is important to note that in these works the initial 
noise was taken to be a classical wave at a specific fre- 
quency and the analysis therefore involves a discrete set 
of modes. 

In the present work we study both theoretically and 
experimentally the growth of the higher order harmonics 
of MI in the regime where the instability arises by spon- 
taneous amplification of vacuum fluctuations. The essen- 
tial difference with respect to the earlier work mentioned 
above is that in the present case a continuous spectrum 
of vacuum fluctuations is involved (rather than a discrete 
set of modes). We will characterize each harmonic by 
its detuning, its intensity and growth rate, and its spec- 
tral width. Several of these features would be different, 
or would simply not appear (in the case of the spectral 
width), if one was dealing with a discrete set of modes. 

Our investigation is carried out in the regime of large 
gain, when many photons have been created in the side- 
bands, but before we reach the regime of pump dcpele- 
tion (ic. before the energy in the side bands becomes 



comparable to the energy in the initial the monochro- 
matic wave). One of the main interests of our work is to 
show how, even though we are working in a semiclassical 
regime where many photons are produced, the fact that 
the instability is seeded by vacuum flucutations leaves a 
distinct signature in the spectrum of the sidebands. 

Our analytical predictions are well reproduced in our 
experimental demonstration, based on light propagation 
through optical fibers. This system is particularly well 
suited for the such investigations because light propaga- 
tion in low loss silica fibers is accurately described by the 
Non Linear Schrodinger Equation (NLSE), hence suitable 
to precise theoretical modelling, and because it is fairly 
easy to work in a regime where initial noise is dominated 
by vacuum fluctuations. In fact the appearance of the 
first harmonic -probably seeded by vacuum fluctuations- 
was already reported in the first experimental investiga- 
tion of MI in optical fibers Thus our approach pro- 
vides quantitative explanation for the complex spectra - 
readily observable in experiments- that arise in MI seeded 
by vacuum fluctuations in the regime of large gain. A 
related comparison, but between numerical simulations 
and experimental results- can be found in |ld| . 

Finally we note that the ideas and methods developed 
here are not restricted to the problem of scalar MI, but 
should apply mutatis mutandis to other forms of MI such 
as vectorial MI, and more generally to any unstable dy- 
namical system involving a continuous set of modes (i.e. 
described by partial differential equations), when the in- 
stability is seeded by vacuum fluctuations. We also note 
that our analysis can also be applied if the initial noise 
is classical white noise. 

Section III Al contains our theoretical analysis. Our 
key result is eq. (|16|) which encapsulates our predictions 
concerning intensities, growth rates, frequencies, spectral 
widths of the harmonics. A discussion of the predictions 
provided by eq. Ijltif) . of the hypotheses that go into its 
derivation, and of its interpretation, are given in section 
III Bl Some background material, and details of some 
calculations have been relegated to Appendices. Section 
describes our experimental setup and shows that the 
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experimental results are in good agreement with the the- 
oretical predictions. 



II. THEORETICAL ANALYSIS 
A. Harmonics of Modulational Instability 

Let A(x, t) be the slowly varying field envelope of a 
light pulse propagating in an isotropic single mode fiber 
with carrier frequency u>n, where x is the coordinate along 
the fiber and t = t — x/v g is the time variable of coor- 
dinates moving at group velocity v g . It obeys the NLSE 
(seeRef. @) 
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where (3 2 is the group velocity dispersion at frequency loq 
and 7 = — 2 i U2UJ A ° — with fin the effective linear refrac- 
tion index for light guided in the fiber, n 2 the non linear 
refraction index, A e ff the effective cross section of the 
fiber. With this normalization |j4(:e,t)| 2 is the instanta- 
neous power flowing through the fiber at position x at 
time t. 

We take the unperturbed solution to be 
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with A a constant and 4>nl = lA^x is the non-linear 
phase. As is well known, when (3 2 < (the anomalous 
dispersion regime), the continuous solution eq. J5J is un- 
stable. We give the "standard" derivation of this instabil- 
ity (following 0) in Appendix IaI where we also include 
the predictions of quantum theory (when the instability 
is seeded by vacuum fluctuations) and discuss the regime 
of large gain. 

In the present work we use a slightly different approach 
introduced by Hasegawa and Brinkman Q. This will 
allow us to derive many of the properties of the harmonics 
of MI (which are inaccessible using the standard method 
of Appendix lA^I . 

Our starting point is to parameterize the solution as: 



A(x, t) = A a e i +"'- (1 + e(x, T )fl 2 e}< x ^ 
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where e and a are real variables. Upon insertion of this 
ansatz into the NLSE one obtains the equations 
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2 V4(l + e) 2 2(1 + e) 

-ci c a + 7 |Ao| 2 e = 0. (4) 

The key to the approach of Hasegawa and Brinkman is 
to keep the non linear dependence of A in e and a in 
cq. ©, but to linearise the eqs. (0J. Obviously this 
does not provide an exact solution to the NLSE, but it 



provides an approximate solution which captures many 
of the features of an exact solution (illustrated for in- 
stance by the good agreement between these theoretical 
predictions and our experimental results). 

Upon linearisation of eqs. (0J one obtains the following 
simplified set of equations: 
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The Fourier transforms of e and a are: 
1 
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When (3 2 < and < lu 2 < A^A\l\f5 2 \ these equations 
possess exponentially growing solutions: 



e(x,uj) = e + (uj)e 9X + e_(w)e 
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where 
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As discussed in Appendix ^ one can also write the 
solution to the NLSE as 



A(x,t) 
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where a u are the positive frequency components of the 
initial noise. The importance of this decomposition is 
that in the quantum theory a w should identified with the 
Heisenberg destruction operators. To obtain the relation 
between e± and a we linearize eq. (j3J to obtain 

A= A e l<l,NL (l + e/2 + ia) 

and then compare the two solutions. The details of this 
comparison is given in Appendix iBl 

In the present work we are interested in the regime of 
large gain when gx » 1. This implies several simplifi- 
cations. First we can neglect e~ gx with respect to e +gx . 
Second we note that g has a maximum at 
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hence we only need its value in the vicinity of uj max 

g(uj) ~ g max - |/3 2 |(o^ - u max ) 2 (10) 

with g max — ^Aq. Finally all other functions can be 
approximated by their value at uj m ax (since they do not 
appear in exponentials, but only as prcfactors). 
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With these simplifications, Eqs.[7|lcad to the following 
relation between tr and e 



and the relation between e± and a, derived in eq. (|B1J| . 

becomes 



(12) 
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The condition that gx » 1 implies that many pho- 
tons are created in the sidebands. Thus we can neglect 
quantum ordering problems and we can take a u to be 
classical white noise with moments : 



^<5(wi - . . . 8(u n - J a ( n) ) (13) 



where the sum is over all permutations a of {1, . . . , n}. 

In summary we have obtained an approximate solution 
of the NLSE 



.4 



A e t4>NL (l + e) 1/2 e le/2 



(14) 



with e given by eq. (|12|) and a u classical white noise as 
described in eq. I|13|) . 

In order to study the power spectrum of the harmonics, 
we expand eq. I|14|) as a series in e to obtain 
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where c„ are the Taylor series coefficients 
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T„ (a) = 1 x a X (a — 1) . . . (a — n+ 1) 



the first few coefficients of which are 

C = 1,01 = ^,02-=^ 



4 i ^3 ^ ) ^4 

Recall that e contains both the frequencies around 
+w rnQ2; and —u> max , see eq. <|12ll . Hence e" con- 
tains frequencies around nuj maxi (n — 2)ti) maK , — (n — 
2)u> m ax, —nwrnax- Thus each power in eq. Ijl5(l gives rise 
to a new harmonic around frequencies ±nw mal . Also 
there is a component around u> — which first arises at 
order n = 2. 

To compute the power spectrum of these harmonics we 
take the Fourier transform at frequency noj ma x + 8 of the 
n'th order term in eq. (|15|) : 

A„(n^ ma:c + 8) = ^J dTe^^+ 5 ^c n A Q e^-e n (r). 
The spectral energy density at frequency nu> m ax + 8 is 

+ S) = (\A + 5)?) • 
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The computation of P n (nuj max + <5) is somewhat tedious 
and details are given in Appendix ]B\ The result of the 
computation is: 



P n (nuj n 
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with 



1,2,3,... and the first few Taylor coefficients 



are |c | 2 = 1, |d| 2 = 1/2, |c 2 | 2 = 1/8, |c 3 | 2 = 1/144, 
|c 4 | 2 = 1/1152. 

The appearance of 8(0) should be interpreted as usual 
in this kind of calculation as 8(0) = T/2ir with T the 
duration of the light pulse. Dividing by T yields the 
spectral power density. 

For n = 1 this describes the growth of the MI side- 
bands; and for n > 1 it describes the growth of the har- 
monics. As mentioned above the first harmonic (n = 2) 
also contains a component around lu = 0. One can repeat 
the above calculation to find that it is equal to 4P2, ie. it 
has exactly the same shape as the first harmonic around 
2uj max except that it is centered around tu = and is 4 
times more intense. 



B. Summary of predictions and interpretation 

Let us first summarize the main predictions contained 
in eq. I|16|l : we will then discuss the interpretation and 
limitations of our theoretical method. These predictions 
are: 

1. We recover the well known result that the funda- 
mental instability appears at frequency ±w max = 
±y/2jAl/\f3 2 \, grows at a rate 2g max = 2^\A \ 2 , 
and has spectral width Auj = u\ = l/2y/\02\x- 

2. The power spectrum exhibits a series of Gaus- 
sian peaks centered on frequencies nu) max , n 

.... —3, —2, —1, 0, 1, 2, 3, ... (n = ±1 corresponds to 
the fundamental instability just mentioned). 

3. The n — l'th harmonic appears at detuning uj n = 
±nu m ax, grows at a rate 2g„ = 2ng rnax = 
2wy\Ao\ 2 , and has spectral width Aw — a n = 
y/nai = y/n/2^\/3 2 \x. 

In addition the first harmonic has a component 
around oj = which has the same properties as the 
component around 2cu max except that its intensity 
is 4 times larger. 

4. Whereas the fundamental instability can be stim- 
ulated by a classical signal around the frequency 
^max, the harmonics cannot be stimulated. They 
are entirely determined by the initial noise around 
the frequency of the fundamental instability. We 
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have verified this prediction experimentally by in- 
jecting a classical signal at frequency 2uj max and 
checking that it does not affect power spectrum of 
the harmonics. 

5. The exact intensity of the harmonics is highly sen- 
sitive to any classical noise present initially around 
the frequency uj max of the MI. Indeed the inten- 
sity of the n — l'th harmonic will be multiplied by 
(1/2 + n c i ass ) n where n c i ass is the number of clas- 
sical noise photons per mode initially present. The 
sensitivity thus increases with the order of the har- 
monic. 

The derivation of eq. I|16[) raises some interesting 
points. First of all, the above analysis can be thought 
of as an expansion valid when there is no pump deple- 
tion, ie. when the power in the modulational instability 
sidebands is small with respect to the power in the pump 
beam. When pump depletion can no longer be neglected 
our computation is no longer valid, since we supposed 
that the pump amplitude is constant. Furthermore in 
the regime where the pump gets depleted there will be a 
back-action of the higher order harmonics on the lower 
order ones, an effect we neglected (for instance we only 
considered the contribution of A n to P n , and neglected 
the contributions of A n i, n' > n). 

The ratio 77 between the energy in the fundamental 
sidebands at frequency io max and —0J m ax, and the pump 
energy can be explicitly computed to be 

J dSPi(uj max + S) 

V = 2 TAJ 

hu e 29 ™* x 
' AA 2 ^\^ 

where the factor of 2 takes into account that there are 2 
sidebands. 

The ratio of the energy in the n— l'th harmonic to the 
energy in the pump beam can then be written as 

2 jd5P n {. x + 5) =2 ^ cn? 
1 A Q 

which shows that up to a slowly varying factor |c„| 2 the 
intensity in the n-l'th harmonic is proportional to rf 1 . In 
our computation we assumed that each successive har- 
monic is smaller than the preceeding ones, ie. P n < P n -i- 
This allowed us to neglect the back-action of higher or- 
der harmonics on lower order ones. Obviously this cor- 
responds to the condition r\ << 1. 

The above analysis does not provide an exact solution 
of the Non Linear Schrodingcr equation Q. Rather we 
have only solved the linearised equations of Hasegawa 
and Brinkman cq. Ipjl. Thus our main prediction eq. 
(|16f) cannot be exact. We expect that the growth rate and 
spectral widths of the harmonics are robust predictions 
because they depend only on the fact that the n — l'th 
harmonic is proportional to e". On the other hand the 



prefactor, and in particular the value of the coefficients 
|c„| 2 cannot be predicted correctly in the present ap- 
proach. (Indeed by using as ansatz a non linear function 
different from eq. J2J), all our results would be unchanged 
except the coefficients c„ which would change). Never- 
theless the simple ansatz cq. © gives surprisingly good 
predictions for the intensities, see the experimental re- 
sults reported in section ILTT1 

Let us conclude by sketching how one could carry out 
a more systematic approach to the harmonics of the MI 
that would predict correctly the prefactor of eq. 1(Pd|| . 
(In fact during our first investigations of this problem 
we adopted this approach, but then switched to the ap- 
proach of Hasegawa and Brinkman which is much simpler 
mathematically). This systematic approach is neverthe- 
less interesting, if only because it gives a different point 
of view to the problem. 

Its starting point is the standard approach of the NLSE 
based on the linear equations described in Appendix 1X1 
There we took an ansatz of the form A = (Ao + Ai)e l ^ NL 
and linearised the equations in A\ to obtain eq. lAlfl . 
But this is only an approximate solution. We can take 
into account systematic corrections to this solution by 
considering the ansatz A = (Aq + A\ + A 2 )e l ^ NL where 
A\ is the solution of the linearised equations. One then 
obtains for A2 the equation: 



id x A 2 = ^-d 2 A 2 - 1 A 2 {A 2 + A* 2 ) - 7 A (2|A 1 | 2 + A\). 

(17) 

This is a linear equation for A 2 with an independent 
term. Thus the solution of the equation for A 2 is a solu- 
tion of the homogeneous equation plus a particular solu- 
tion of the inhomogeneous equation. The independent 
term is quadratic in A\ and therefore is proportional 
to e 29maxX and contains frequencies around 2uj max and 
around u> = 0. Upon solving for A 2 one will find that 
A 2 is largest around 2ui max and around lo = and is 
proportional to e 29mamX . Thus A 2 will encode the be- 
havior of the second and zero'th harmonic. Successive 
orders in perturbation theory will give rise to the suc- 
cessive harmonics. Wc leave the detailed investigation of 
this approach to future work. 



III. EXPERIMENTAL RESULTS 

Our experimental set-up is reported in Fig. ^ It con- 
sists of a Q switched laser (Cobolt Tango) that produces 
pulses at 1536 nm, with a 3.55 ns full-width-at-half- 
maximum duration T and a 2.5 kHz repetition rate /. 
The pump power Pq = Aq is adjusted using a variable 
neutral density filter (ND). A narrow band filter (NBF) 
allows a wide spectral range around the pump wavelength 
to be free of noise. A polarizing beam splitter (PBS) en- 
sures that the pump is linearly polarized and allows to 
measure the injected power proportional to the rejected 
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FIG. 1: Experimental setup. ISO: isolator, ND: variable neu- 
tral density filter, NBF: narrow band filter, PBS: polarizing 
beam splitter, -|: half- wave plate, PMF: polarization main- 
taining fiber, PBG: fiber Bragg grating, and OSA: optical 
spectrum analyzer. 
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beam. We used the Fibercore HB1250P polarizing main- 
taining fiber, although the experiment could equally have 
been realized with non birefringcnt fiber. A half-wave 
plate is used to ensure that the pump polarization is 
aligned with a principal axis of the fiber, whereupon po- 
larization effects can be neglected and the above results 
for scalar MI apply. The fiber parameters are deduced 
from both scalar and vector modulation instabilities (see 
Ref. 0])- The fiber length L is 51 m. The group- velocity 
dispersion parameter /?2 is —15.27 ps 2 km -1 . The Kerr 
nonlincarity parameter 7 is 3.26 W~ 1 km~ 1 . (The beat 
length, which is irrelevant to the present experiment, is 
17.9 mm). Lastly, a fiber Bragg grating (FBG) rejects the 
pump wavelength before the measurement of the spectra. 
The rejection of the pump avoids detector blinding and 
allows us to reach the sensibility limit of the optical spec- 
tral analyzer (OSA). 

The Narrow Band Filter (NBF) eliminates all resid- 
ual photons except those at the pump wavelength. This 
ensures that the MI and harmonics indeed arise from 
vacuum fluctuations and not from classical noise. We 
have checked that this is indeed the case in two ways. 
First numerical simulations of the Stochastic Non Linear 
Schrodingcr equation reproduce very well the observed 
spectra, including the harmonics [Hj. Second we have 
studied in detail, using a single photon detector, the 
MI in the regime where relatively few photons are pro- 
duced. This investigation shows that there is in fact a 
small amount of noise present due to spontaneous Ra- 
man scattering in the fiber. But the number of Raman 
photons per mode is much smaller than 1 which means 
that they do not affect the spectra when the gain is large. 
In summary both investigations show that the MI pro- 
cess is dominated by spontaneous effects and stimulation 
by classical noise is negligible. 

A sample of the collected spectra is shown in Fig. |3 
for different pump powers Pn- We clearly observe the 
growth of the MI and the appearance and growth of the 
harmonics as the pump power increases. The largest pair 
of peaks correspond to the MI. Then each new harmonic 
gives rise to a pair of peaks further and further from the 
pump wavelength. In addition there is a peak around 



FIG. 2: Spectra of the output field for increasing peak power 
P equal to 31.6 W, 38.9 W, 45.9 W, 49.1 W, 52.2 W, 57.4 W, 
59.6 W, and 62.4 W. The resolution bandwidth is 0.1 nm. 
Note that the Fiber Bragg Grating has removed the pump 
wavelength, thereby avoiding detector blinding and allowing 
us to reach the sensitivity limit of the OSA. 



the pump wavelength. Because this peak overlaps with 
the pump it is largely rejected by the FBG and cannot 
be well characterized from this first set of measurements. 
In consequence we have characterized the peak around 
the pump wavelength from a second set of measurements 
obtained without the FBG (figure not shown). For each 
pump power Pq , each harmonic has been fitted according 

to a gaussian function 7 n (Pn)exp ^— ^"^"(Pb)^ ) wnere 
/„, Lu n , <j„ characterize the intensity, frequency and spec- 
tral width of the harmonic. 

The measured intensities I u (Pq) of each harmonic are 
compared in Fig. to the theoretical predictions (con- 
tinuous curve). To obtain the theoretical curves, Eg. 1161 
must be integrated over the OSA resolution bandwidth 
(0.1 nm). Moreover we must identify 5(0) = T 2 * p^/ , 
where T is the pulse full width at half maximum, ST 
is the sweep time equal to 27.9 s, / is the laser repeti- 
tion rate, and PTS denotes the number of samples per 
spectrum equal to 1001. In order to obtain a good over- 
lap between theoretical curves and measured points, we 
had to further shift 8(0) by 7 dB. This discrepancy is 
not unreasonable given that the above theory was based 
on a continuous pump, whereas now we are dealing with 
a gaussian pulse, and given our inability to carry out 
absolute measurements with the OSA. With this shift 
the theoretical and experimental intensities agree well. 
This shows that eq. I|16|) , including both the exponential 
terms and the prefactor, correctly predicts the absolute 
intensities of the harmonics. We further confirmed this 
agreement by fitting the intensity of each harmonic to an 
exponential law I n (Po) ~ A n exp(a n 2"f P L) where A n is 
an n dependent constant, and a n a growth factor (dotted 
lines in Fig. 0). Theory predicts that a n /ai = n. This 
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FIG. 3: Maximum intensity as a function of the peak power 
of the pump. The up-triangles, down-triangles, circles, and 
squares represent respectively the fundamental MI peak (n = 
1), the first (n = 2), second (n = 3) and third (n = 4) har- 
monics. The dashed lines are a fit to the data according to 
an exponential law. The continuous curves are the theoretical 
predictions of Eg. 1161 As described in the text the experimen- 
tal points were shifted by 7dB to agree with the theoretical 
predictions. For clarity, we show a spectrum. 



prediction is well verified, see Table Q] 

Note that around 60 W the total energy in the side- 
bands becomes comparable with the pump energy: the 
MI saturates and our analysis is no longer valid. 

Theory predicts that the harmonic around ui = 
should have the same growth rate, but 4 times the in- 
tensity, as the harmonic around 2u max . There is reason- 
able agreement concerning the growth rate, see Table [I] 
We were unable to check the factor of 4 as the intensity 
measurements of the harmonic around u> = were not 
precise enough, although it is always more intense than 
the harmonic around 2io max , see Fig. [21 

We now turn to the angular frequency shifts lo„(Pq). 
In Fig. 21 we plot the measured values and the theoreti- 
cal predictions (continuous curves). We note a very good 
overlap. Moreover, from the data in Fig. 01 we have com- 
puted the ratios between the (n— l) th harmonic frequency 
Lu n (Po) and the fundamental frequency u>i(Pq). The av- 
erage values (u) n /u>i)p g (where the average is over the dif- 
ferent values of Po) are reported in Table|U There is good 
agreement with the theoretical prediction u) n /uj\ = n. 

Finally, the measured spectral widths a n are reported 
in Fig. [S] as a function of Pp. T hey are compared to 
the theoretical values y / n/4|/32|i (straight lines). The 
agreement is good for the fundamental frequency (n = 1) 
and the first harmonic (n = 2). In the case of the second 
(n = 3) and third (n = 4) harmonics it is not possible 
to conclude, although the values are consistent with an 
increase in spectral width with the order of the harmonic. 
This is because for high pump powers (> 60 14^) the MI 
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FIG. 4: Angular frequency shift as a function of the peak 
power Po of the pump. The up-triangles, down-triangles, cir- 
cles, and squares represent respectively fundamental MI peak 
(n — 1), the first (n = 2), second (n = 3) and third (n = 4) 
harmonics. The continuous curves correspond to the theoret- 
ical prediction uj n = n-\/2-fPo/\ft2\- 
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FIG. 5: The spectral width a„ of each harmonic is shown 
as a function of the peak intensity Po. The up-triangles, 
down-triangles, circles, and squares represent respectively the 
fundamental MI peak, the first, second and third harmon- 
ics. Straight lines correspond to the theoretical prediction: 
a n = y/n/4\p 2 \L. 



saturates; in this regime we find that the spectral width 
increases for all values of n. And when the peaks of 
the MI harmonics have low intensity they are broadened 
for an unknown reason. For the second (n = 3) and 
third (n = 4) harmonic the widths are always affected 
by either of these effects. Finally the ratios (cn/ a i)p 
(where the average is taken over those values of Po where 
a n is approximately constant) are reported in Tab.[I] As 
expected, only the first harmonic leads to a proper value. 



7 



TABLE I: Behavior of the harmonics. * and f denote respec- 
tively the first and the second set of measures. The theoretical 
prediction is that all the quantities are equal to n (except the 
center frequency of the second harmonic centered on w = 0). 



n (Order of 
harmonic = n — 1) 


Measured values 




2 (w ~ 0) 
2 (lo ~ 2u max ) 


2.11±0.03 t f 2.43±0.15 f 
2.31±0.02 t 1.95 ±0.01+ 1.99 ±0.12+ 


2 (u ~ 2a; maa! ) 
3 
4 


2.15 ±0.07* 1.85 ±0.01* 1.77 ±0.06* 
2.21 ±0.06* 2.95 ±0.01* 2.92 ± 0.23* 
4.89 ±0.15* 3.94 ±0.03* 3.36 ± 0.61* 



IV. CONCLUSION 

In summary we have studied the higher order harmon- 
ics of Modulational Instability in the regime where the 
MI arises spontaneously through amplification of vacuum 
fluctuations. We have shown that in this regime there is 



a rich phenomenology which can be predicted theoreti- 
cally and is well verified experimentally In particular we 
have obtained predictions for the dctunings, intensities, 
growth rates and spectral widths of the harmonics, see 
our key result eq. (|16f) : and found good agreement with 
our experimental results based on light propagation in 
optical fibers. We note that the shape of the spectrum, 
and in particular the relative intensities of the higher or- 
der harmonics, is highly sensitive to the initial presence 
of classical noise, and can therefore be taken as a signa- 
ture that the MI is seeded by vacuum fluctuations. It 
should be possible to extend our work in a number of 
directions, for instance to predict exactly the coefficients 
\c n \ 2 appearing in cq. I|16(l . to study the regime where 
the pump starts to get depleted, and to extend it to other 
forms of MI, or to other kinds of instabilities. 
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1. Linear Perturbation Theory 

In order to solve the NLSE eq. (|TJ we consider an 
ansatz of the form 

A= (A Q + A 1 )e l ' t ' NL . 

where A\ is a small perturbation. We linearise the NLSE 
around the continous solution to obtain the equation for 
Ai. 



id le A 1 = ^%A 1 -<yA%(A 1 +Al) . 

We then carry out a Fourier expansion of A\. 
1 



Ai{t,x) 



>2tt 



du> ax{x, uj)e 



(Al) 



(A2) 



where cu is the detuning with respect to ujq. Eq. (|A1|) 
then becomes 



APPENDIX A: MODULATIONAL INSTABILITY 



In this Appendix we recall the usual approach to MI. 
We first solve the linearised equations, then consider the 
predictions of quantum theory, and finally consider the 
regime of large gain. This appendix thus contains back- 
ground material useful for reading the rest of the article. 



id x a\(x,uj) = ~^-aj 2 ai(x,uj)-jAl (ax(x,u) ± a*(x, —u)) 

(A3) 

When /3 2 < and when < u 2 < uj 2 = 4 7 ^/|/3 2 | 
this equation possesses exponentially growing and expo- 
nentially decreasing solutions. On the other hand when 
u! 2 > lo 2 or when fa > cq. (|A3(I possesses oscillating 
solutions. Let us consider the unstable solutions. These 
are proportional to e ±gx where the gain is 



2 I 



1/2 



(A4) 
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We can reexpress the solution in terms of the initial 
conditions as 



a\{x,oj) = ai(0,u>)/j,(x,u)) + a*(0, w) (A5) 



H(x,w) = cosh(gx) + i ^ A2 ° + sinh(gx)(A6) 



v(x, u>) 



fj 



■ six\h(gx), 



(A7) 



where fi and v obey the condition \/j,\ 2 — \v\ 2 = 1. 



2. Quantum Theory 

The above solution of the NLSE can also be used to 
derive the form of the instability when the Modulational 
Instability is seeded by vacuum fluctuations. In this case 
cq. should be interpreted as the Heisenberg equa- 
tion for the operators A and A* . The Fourier transform 
cq. (|A2() defines the destruction operators ai(x,ui). Their 
hermitian conjugate a\(x,u>) are the creation operators. 
These operators obey the commutation relations 

[Si (x,uj), a\(x, u')] = h(u>o + oj)S(uj — u>') . (A8) 

Equation (|A5|) should then be reinterpreted as giving the 
relation between the creation and destruction operators 
at distance x along the fiber and at the origin x = 0: 

ai(x, w) = Si(0, co)/j,(x, ui) + a\(0, w)v(x, u>) . (A9) 

The initial state |0) contains no photons: 

ai(0,w)|0) = . (A10) 



From eq. (|A9|) one can then compute the expectation 
values of products of creation and destruction operators 
at distance x along the fiber. For instance one finds 

(0\a\(x,uj)ai(x,uj')\0) = Ti(ujq — ui)\i/(x, lo)\ 2 8(lo — u>'). 

(All) 

From this it follows that the spectral power density at 
frequeny u>q + to, P(x, u>), is given by 

P(x,u) = (0|a t 1 (a;,w)ai(.T,w)|0) 

= h{uj a - lu)\v(x,lu)\ 2 5(0) 

~ huj a \is(x,u})\ 2 S(0). (A12) 

It is infinite because we have computed the power when 
the pump beam is monochromatic and of infinite dura- 
tion. If we suppose that the pump beam lasts only for 
a duration T then one should -as usual in these kinds of 
situations- interpret this infinity as 



T 

*<B>- S . 



(A13) 



Then the spectral power density at position x at fre- 
quency Do + w is P(x,u>) = Tiloq\v\ 2 T/2'k. 

Note that the quantum solution also predicts other ef- 
fects such as correlations (two mode squeezing) which 
will not be studied here. 



3. Behavior in the regime of large gain 

We now consider the regime where the gain e gx is large. 
In this case we can simplify the solutions obtained above. 
First of all we note that the produced photons will be 
localised around the frequency uj max at which the gain is 
maximum: 



2l3 
\(h\ ^ 



The maximum gain is given by 



9n 



■ g(iOmax) = "fAl 



(A14) 



(A15) 



In order to describe the behavior in the vicinity of the 
maximum gain we can expand the gain around u> max as 



g(u>) ~ g m ax 



y ( x2 

— \ LO — LOmax ) , 



(A16) 



v{x, u) 



I", 

■exp[xg r , 



xg"/2(u;-u ) 2 }, 
xg"/2(u;-u )% 



where g" = 2\f3z\. In the vicinity of uj m ax the coefficients 
fi and v simplify. They can be approximated by: 

1 

2 

i 

2 

where we have kept only the exponentially growing terms 
in v, used the approximate expressions derived in eq. 
(| A 16(1 and dropped the u> dependence of the prefactors 
in /x, v. 

Furthermore in the regime of large gain the quantum 
solution simplifies. Indeed since there are many photons 
in each mode one can carry out a semiclassical treatment 
in which one neglects ordering problems. Thus in this 
regime one can reproduce the predictions of the quantum 
solution by taking the initial conditions of the classical 
solution a(0,<x>) to be white noise with power Tiujq/2 per 
mode. More precisly the a(0,u>) should be taken to be 
complex delta correlated gaussian random variable dis- 
tributed according to the probability distribution 

P(a(0,w)) = J-e-N 2 /^ (A17) 



2-no 2 



with variance a 2 such that 



(a(0,w)a*(0,u/)> = -y-S(u-o/). 



(A18) 



This will correctly reproduce the quantum predictions up 
to corrections proportional to exp[— gx]. Indeed taking 
the probability distribution to be gaussian correctly re- 
produces in the regime of large gain the combinatorial 
factors which arise from Wick contractions when expec- 
tation values of the products of many creation and de- 
struction operators are taken. Thus for instance in this 
regime we have 

(o(0, wi) . . . o(0, uj n )a* (0, . . . o*(0, <)) = 

9 » 8 n ,m 6 ^ - ^(1)) • • • 5 ^ n ~ W ^(n))' ( A19 ) 
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where the sum is carried out over all permutations a of where for conciseness we note = ai(0, w). This term 



In the regime of large gain the spectral power density 
is thus 

P(x,uj) = (a* (x,uj)a(x 1 uj)) 
T Tilo 

- cxp[25f ma3; x] exp[-#"(w - u max ) 2 x] . 



2tt 4 



(A20) 

APPENDIX B: DETAILS OF CALCULATIONS 

1. Relation Between e, a and a 

Upon linearising eq. (j3J) in e and a we obtain 

A = A e ilt > NL {l + 6/2 + icr) 

Upons inserting the forms given by eqs. (jSJ) and Q one 
finds 



contains frequencies around: 

ncj max , {n 2)u^ maa; , . . . , {ri 2>^uj rnaX: Jicj max . 
The only new frequencies which appear at order n are 
thus around ±nu) max (except for n = 2 when new fre- 
quencies appear also around u! = 0). It is the behavior 
of the n'th order term around these frequencies which we 
are interested in. 

We thus take the Fourier component of c n e tl >' NL Aqc 71 
at frequency nuj max + 6 to obtain 



A n (x,nLu max + S) 



CnApe^L 
( 27r )(n-i)/2 



n 



- UJm - 3=1 

«w mai - S 



A = A e l 



1 



1 



27T Jo 



1 



\ j=1 
x exp 

n 

><n 

3=1 



X^( w 3> 

V3 = l 



An 2 



A, 



'5_ 

'5 



£+ H(- + — ^)e— +^ 
2 |/? 2 |t^ 



(B3) 



1 

- l A 2 \I3 2 \UJ 2) > 



The spectral density of power at frequency nuj max + 5 is 
P n (nuj max + 5) = (\A n (nu max + S)\ 2 ) 
\A \ 2 \c n \ 2 (RtJo)" 



(27T)"- 1 2"|A | 2 " 



3=1 3=1 



We can then identify 



S io 3 ■ - nu max - 5 \ 6 w'j 



noj Tl 



^3=1 



^'=1 



e+(w) 



e_M = 



ai(0,w) / ,|/3 2 |^ 2 \ , aJCO.w) / , .|/3 2 |w 



2An 



1 - i± 



2g 



2A () 



l + i 1 



'Mi 



exp +ffK) a 



oi(0,o;) / , .|/3 2 |w 2 \ al(0,a;) / 



2A 



1 + i- 



2<? 



2A, 



1 - i- 



2.9 



v i=i 



(Bl) 



— £ 2 x - w ; (1) ) ... 2 x 8(u n - w ; (n) ) ,(B4) 



which in the regime of large gain reduces to eq. I|12|) . 



where the sum over a is a sum over all n\ possible permu- 
tations. Carrying out the integrals over ui'j this becomes: 



2. Power Spectrum of the Harmonics 

Here we give the details of the calculations leading from 
eq. (|15|) to eq. (|16f> . If we explicitise eq. I|15|) we obtain: 

A n (x,r) = c n A e l ^ L e n 
c n A e^t- 
(2tt)"/ 2 

~[ dujj exp -i ^2 exp ^ g{uj)z 



P„(nw TOax + <5) = 5(0) 



|Ao| 2 |c„| 2 (fia; )"n! 
(27T)"- 1 2"|A | 2n 



3 = 1 



ME 



cj, — mo, 



J =1 



tax 

S exp 2^2g(ujj)x . (B5) 

3 = 1 



To evaluate the remaining integral we make the change 
of variables 



3 = 1 
n 

n 

3=1 



3=1 

a-uj, 1 — i a *-uj, l + i 



,3 = 1 



LOj — UJmax H h £ 

n 

In terms of these variables we have 



An 2 



An 



(B2) 



3=1 



3=1 
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(since X)?=i £7 = 0)- We recall that g" is equal to 2|/?2|- 
We thus have 



/n / n 



Wi — TO„ 



<5 ] exp [ 2^g(wj)a; 



exp (2ng ma:c a:) exp(-g"5 2 a;/n) 



n-l 



cxp[2ng maa; x] exp[-g"xS /n] (B6) 



where the integral over £j is carried out as follows: change 
variables to (j = Rjk^k where R is an orthogonal 
matrix such that £i = X)fc Cfc/V^- Then the Jacobian of 
this transformation is 1. The integrals over Q factorise 
into one delta function and n — 1 gaussians, yielding eq. 
(|B6|) . Inserting this into eq. (|B5(1 yields eq. JT^Jl. 



